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ABSTRACT 


An  improved,  nonlinear,  constrained  mathematical 
programming  optimization  algorithm  is  presented  in  this 
report.  It  couples  a rotating  coordinate  pattern  search 
with  a feasible  direction  finding  procedure  used  at  points 
of  pattern  search  termination.  The  procedure  is  compared 
with  nineteen  algorithms,  including  most  of  the  popular 
methods,  on  ten  test  problems.  These  problems  are  such 
that  the  majority  of  codes  failed  to  solve  more  than  half 
of  them.  The  new  method  proved  superior  to  all  others  in 
the  overall  generality  and  efficiency  rating,  being  the 
only  one  solving  all  problems.  It  was  particularly  effec- 
tive on  constrained  problems  where  it  was  best  in  all 
rating  categories. 
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NOMENCLATURE 


behavior  of  constraints  g^(x^) 

e^  “ small  arbitrary  constants  defining 
finite  difference  step  size  used  to 
calculate  the  gradient  of  the  objective 
function,  nearness  to  constraints,  and 
convergence  criterion  respectively. 

efficiency  ratings  (eqs.  30  and  31) 
objective  function 

composite  objective  function 
constraint  function 

number  of  variables 
number  of  constraints 

set  of  active  lower  and  upper  regional, 
constraints,  respectively 

lower  limit  on  the  behavior  of 

generality  ratings  (eq.  29) 

penalty  function  used  when  a constraint 
Is  violated 

best  feasible  direction  vector 
generality  and  efficiency  rating  (eq.  32) 
upper  limit  on  the  behavior  of 
problem  variables  or  variable  vector 
optimal  variable  values 

weighting  factor  for  constraint  g^(x^) 

used  in  the  feasible  direction  finding 
problem 


a 


step  size 
minimum  step  size 


■ constraint  activity  limits  for  the  direction 
finding  problem 

■ objective  function  and  variable  in  the 
feasible  direction  finding  problem 

■ gradient  operator 


magnitude  of  <j> 


Superscripts 


lower  limit 


* at  rth  iteration 


- local  exploration  base 

■ transpose  of  vector 
* upper  ,limit 

■ comparison  quantity 
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INTRODUCTION 


A variety  of  methods  involving  ordinary  and  variational 
calculus,  mathematical  programming,  and  a number  of  special 
techniques,  such  as  the  fully  stressed  concept  used  in  struc- 
tural design,  are  available  to  treat  optimal  design  problems. 
Among  these  methods,  the  mathematical  programming  procedures 
appear  to  have  the  broadest  range  of  application  [1]^.  Such 
methods  are  flexible,  easy  to  adapt,  and  can  offer  "automatic" 
optimal  computer  solutions  [2], 

Mathematical  programming  methods  solve,  or  approach  the 
solution  to,  the  problem:  Find  those  values  of  the  variables, 

x^,  that  minimize  (or  maximize)  the  objective  function 

f(x±)  i * 1,2,  ...  I (1) 

such  that  all  constraints 

gj  <xi)  i 0 j - 1,2,  ...  J (2) 

are  satisfied.  The  objective  function  (often  called  the  merit 
or  payoff  function)  quantitatively  defines  the  meiit  of  the 
design  as  a function  of  the  design  variables,  x^.  Inequal- 
ity constraints  are  used  to  control  or  define  the  behavior  of 
the  design  (behavior  constraints). 

Inequality  constraints  can  also  be  used  to  limit  the 
values  of  the  variables  within  specified  limits  (regional  con- 
straints) . 


^"Numbers  in  brackets  designate  References  at  end  of  report. 
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It  is  usually  convenient  to  treat  such  constraints  some- 
what differently  than  behavior  constraints.  Regional  con- 
straints can  be  given  In  the  form 

1 1 (3) 

where  x.^  and  xi  are  the  lower  and  upper  limits  respectively  on 
the  variable  x^. 

A number  of  effective  methods  exist  for  special  forms 
of  this  problem,  such  as  the  unconstrained  problem,  the  linear 
programming  problem,  and  many  quadratic  programming  problems. 
Relatively  few  effective  methods  are  available,  however,  to 
treat  the  form  most  often  encountered  in  design,  the  nonlinear 
constrained  optimization  problem.  Even  the  better  nonlinear 
prgramming  methods  have  limitations  on  the  size  and  complex- 
ity of  problems  they  can  solve  within  a justifiable  amount  of 
computer  time.  Large  order,  highly  nonlinear  systems  requir- 
ing lengthy  design  computations  can  be  quite  formidable.  A 
typical  nonlinear  mathematical  programming  procedure  requires 
several  hundred  to  several  thousand  sets  of  objective  and 
constraint  function  evaluations  to  approach  the  optimum  on 
most  multivariable  problems.  Furthermore,  reliability  Ls  a 
major  problem  with  these  methods  [3] . 

The  new  Direct  Search-Feasible  Direction  (DSFD)  algo- 
rithm appears  to  offer  superior  performance  with  respect  to 
speed  and  reliability  [4],  It  is  the  only  procedure  studied 
that  solved  all  of  Eason  and  Fenton's  test  problems  [3-4]. 
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THE  OPTIMIZATION  ALGORITHM 
The  basic  search  , , 

The  "rotating  coordinate"  (RC)  pattern  search  used  in 
[5]  and  described  in  detail  in  [6]  is  employed  as  the  basic 
optimal  search  technique.  The  usual  pattern  search  strategy 
[7]  is  utilized  in  this  procedure  except  that  after  establish- 
ing a pattern  move,  the  first  local  exploration  step  is  taken 
in  the  direction  of  the  move  and  all  other  local  steps  are 
taken  normal  to  this  direction  rather  than  in  the  direction  of 
the  variable  coordinates,  as  in  the  usual  pattern  search  pro- 
cedure. The  movement  of  this  search  on  a hypothetical  composite 
objective  surface  using  the  penalty  function  of  [5]  is  shown 
in  Fig.  1.  It  may  be  seen  from  Path  1 that  in  the  feasible 
region  the  search  tends  to  move  essentially  in  the  direction 
of  the  gradient,  while  near  the  feasible-infeasible  boundary 
it  moves  easily  along  the  ridge.  The  local  steps  taken 
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normal  to  the  pattern  move  provide  the  turning  component  nec- 
essary to  maintain  efficient  movement.  It  should  be  noted 
that  the  method  does  not  require  a feasible  starting  point  as 
do  most  procedures  (see  Path  2).  The  closed  dots  In  th.is  fig- 
ure indicate  unsuccessful  steps  while  the  open  dots  are  suc- 
cessful moves.  The  digits,  associated  with  some  open  dots, 
identify  the  pattern  base  numbers. 

This  basic  search  procedure  was  found  to  be  substan- 
tially superior  to  the  original  Hooke  and  Jeeves  search  in 
earlier  comparison  studies  [5-7]. 

Since  the  basic  search  scheme  is  formulated  to  treat 
unconstrained  optimization  problems,  the  constrained  problem 
of  (1)  - (3)  is  transformed  to  the  form;  Find 

I 

F^)  - min  F(xi>  (4) 


where 


7(.xt)  * f (x^  + P(x£)  . 


(5) 


Here  P(x^)  is  the  largest  of  the  penalties  Pj(x^)  which 
are  of  the  form 


Pj(xi)  “ “ Xj 

where 

[*]  - * <j>  < 0 

[4>]  - 0 4>  > 0 . 

When  0>  g.  > - e»  then 
J ^ * 

Xj  - ZlfCx^  - f(x±  + Axi)]/[gj  (x±)  - gj(xt  + Ax^I 


(6) 


(7) 


(8) 


! 
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where 


Ax^  - eiVf(x£)/|Vf(x£) |.  (9) 

The  quantity  e^  is  an  arbitrary  small  real  number  representing 

the  size  of  the  step.  Ax.  taken  in  the  direction  of  the  gradi- 

-L  » 

ent  of  fCx^  where  this  gradient  is  evaluated  at  the  local  ex- 
ploration base  point  x^.  it  may  be-seen  that  a penalty  is  used 
only  if  a constraint  is  violated.  If  the  violation  of  any  of 
these  constraints  is  greater  than  an  arbitrarily  specified 
amount  e^,  then  X^  is  not  calculated  from  eqn  (8)  but  is  made 
arbitrarily  large.  That  is,  when 


gj (xi>  < “ e2 


xj  ■ K 


where  K is  an  arbitrary  large  number. 

The  justification  for  the  use  of  this  penalty  function 
form  is  discussed  in  detail  in  Ref.  [5]. 

Violations  of  eqns  (3)  are  handled  by  prohibiting  moves 
that  violate  any  of  these  equations. 

All  of  the  constraints  are  evaluated  only  at  each  tem- 
porary base  (points  reached  after  a pattern  move) . For  all 
other  points  in  the  search  only  nearby  constraints,  this  is 
those  where 

gj  > e2  (H> 

are  evaluated  since  it  is  presumed  that  only  these  constraints 
need  be  considered. 


6 


The  Secondary  Search 

Unfortunately,  the  pattern  search,  even  with  the  rotating 
coordinate  improvement,  is  not  particularly  reliable  [5,6]. 

Thus,  some  method  is  needed  to  determine  if  the  point  of  pattern 
search  termination  is  optimal  and  if  it  is  not,  to  determine  a 
direction  in  which  to  restart  the  search.  Pappas  and  Amba-Rao 
[8]  and  Siddal  [9]  apply  a secondary  search  at  points  of  pat- 
tern search  termination.  Unfortunately,  the  earlier  secondary 
search  strategies  although  more  powerful  locally  than  the  pat- 
tern search  do  not  provide  a rigorous  procedure  for  confirm- 
ing optimality. 

The  direction  finding  procedure  of  Zoutendijk  [10]  can, 
however,  be  used  to  confirm  the  optimality  of  a point  of  pat- 
tern search  termination  and  determine  a direction  in  which  to 
restart  the  pattern  search  if  the  point  is  not  optimal.  The 
direction  finding  problem  can  be  stated:  Given  ttve  set  x^, 

find  the  set  s^  that  results  in  a 


for  which 


max  0 


0 > 0 


(si)iVf(xi)  + 0 < 0 


(Si}  Vgj (xi>  + Wja  1 0 J 


-1  < s.  < 0 k e K“ 
— k — a 

0 < 3k  < ! k e Ka+. 
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Here  (s^)  indicates  the  transpose  of  vector  s^f  is  a 
weighing  parameter,  the  set  J contains  the  active  constraints 

Si 


gj  > " e. 


(18) 


where  e ^ , an  array  of  small  arbitrary  positive  constants  defin- 
ing "activity,"  and  K K+  constitute  the  active  upper  and 
lower  side  constraints,  respectively. 

Zoutendijk  [10]  shows  that  if  the  solution  s.^  is  a null 
vector,  then  the  point  is  a local  optimum  and  if  not  then  the 
direction  of  s^  is  the  best  feasible  direction. 

Equations  (12)-(17)  constitute  a linear  programming  prob- 
lem with  the  variable  Sj,  and  O.  Equation  (12)  is  the  objective 
function  and  the  remaining  equations  the  constraints.  The  solu- 
tion can  be  obtained  reliably  and  efficiently  using  any  suit- 
able linear  programming  method  such  as  the  simplex  procedure. 

The  Optimization  Procedure. 

The  general  procedure  here  is  similar  to  that  used  in 

[3]  and  [6]  except  the  pattern  search  and  optimality  checks 

used  in  [5]  are  replaced  by  those  described  previously.  An 

0 * 

arbitrary  initial  starting  point  x^  is  selected,  F defined 
as  equal  to  F(x^),  and  a quantity  D set  equal  to  unity.  The 
optimal  search  is  then  started  using  the  RC  procedure  with  a 
step  size  and  continued  until  it  terminates  at  some  point 
x^,  where  the  secondary  search  is  invoked. 

If  the  pattern  search  terminates  at  a point  in  the  fea- 
sible region,  that  is,  where  equation  (3)  is  satisfied  for  all 
j,  the  direction  finding  problem  is  formulated  with  = 1 and 
solved  at  the  point  x*.  The  components  of  the  gradients  are 


generated  using  simple  forward  differences  and  thus  no  analyti- 
cal derivatives  are  employed. 

If  the  pattern  search  terminates  in  the  infeasible  region 
at  a point  near  the  feasible-infeasible  boundary,  that  is,  if  some 

e2  - (*£)  — 0 C19) 

(see  Base  15  of  Path  1 in  Fig.  1)  a weighting  parameter  of 

Wj  ■ 100  is  used  to  help  drive  the  search  to  the  feasible  region. 

In  the  case  where  the  pattern  search  terminates  in  the  in- 
feasible region  away  from  the  boundary,  that  is,  where  any  con- 
straint equation  violates  the  left  side  of  (19),  the  search  is 
abandoned  since  in  this  instance  the  pattern  search  has  failed  to 
generate  a feasible  point,  even  with  the  use  of  a large  penalty 
parameter . 

Once  the  direction  vector  is  computed,  the  quantity 
F(x^+1)  is  evaluated,  where 

x*+1  - x*  + Ax±  (20) 

and 

Ax±  - ctrsi/ 1 [ l sm2 ] 1/2  { if  s1  * 0 (21) 

m=l 

Here  ar  is  th.e  step  size  used  for  tlxe  pattern  search  just  termi- 
nated. If  any  constraint  not  considered  in  the  formulation  of 
the  direction  finding  problem  has  been  violated  at  x^  , then 
the  direction  finding  problem  is  reformulated  considering  these 
constraints  and  a new  s^  and  F(x^  ) computed. 

If  8^  « 0 the  activity  definition  limits  and  the  step 
size  ar  are  halved.  If  the  set  changes  as  a result  of  the 
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change  in  due  to  one  of  the  constraints  becoming  deactivated 
as  a result  of  the  reduction  in  the  activity  specification  the 
direction  finding  problem  is  reformulated  and  solved.  Otherwise 
the  activity  limit  and  step  size  is  further  reduced.  This  proc- 
ess is  repeated  until  a nonzero  solution  of  s^  is  obtained  or 


until 


where  a 


ar  < a . 

— min 


(22) 


is  a specified  minimum  basic  step  size. 
Now  if 

F(*^+1)  < F<xJ) 

. r+1 


(23) 


the  RC  search  is  restarted  using  x^"r'L  as  th_e  new  base  point 
in  the  pattern  search  strategy.  If  equation  (22)  is  not  satisfied, 
then  a is  refined  as 

<xr+1  ■ or/2 


Now  if 


or  if 


„r+l  „ 
a <ot 


min 


(24) 


(25) 


and 

where 


»r  < «3 

Dr<  D*  ^ (26) 

Dr  - [F*  - F(x*)  ]/F* 

and  where  «min  and  e^  are  arbitrary  small  positive  variable  and 
objective  function  convergence  constants,  respectively,  the  search 
is  abandoned  and  the  point  is  assumed  to  be  sufficiently  near  opti- 


mum. Otherwise  a new  x 


r+1 


is  defined  by  equations  (20-21)  with 


or  replaced  by  ar+^  until  equation  (23),  (25),  or  (26)  is  satisfied, 
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Once  equation  (23)  is  satisfied,  then  the  RC  search  is  restarted, 
with  x*+*  as  a new  base  point  in  the  pattern  search  strategy, 

* * ft  j* 

with  step  size  a , F ■ F(x£  ) , and  D - D . If,  however,  two 
successive  step  size  reductions  fail  to  satisfy  equation  (23), 
then  the  entire  optimization  procedure  is  restarting  from  point 
x*»  with  step  size  or+^>F  ■ FCx*)^  and  D ■ Dr  and  repeated 
until  equation  (25)  or  (26)  is  satisfied. 

It  may  be  seen  that  the  step  defined  by  equation  (21), 

% 

taken  when  s^  is  nonzero,  is  a move  in  the  best  feasible  direc- 
tion. Thus,  if  equation  (23)  is  not  satisfied  by  this  move,  it 
indicates  that  the  step  size  is  too  large.  The  step  size  is  then 
reduced  in  an  effort  to  locate  a better  point,  first  by  taking  a 
smaller  step  in  the  s^  direction,  and  in  the  event  this  fails, 
by  attempting  to  restart  the  pattern  using  this  smaller  step. 

The  process  is  repeated  until  convergence  of  the  objective  func- 

a 

tion  is  achieved  or  a minimum  step  size  is  reached. 

A zero  s vector  indicates  the  presence  of  an  optimum. 
Since,  however,  in  the  consideration  of  the  active  constraints, 

Ej  V 0>  the  point  x*  may  be  merely  near,  rather  than  at,  the 
optimum.  Thus,  the  activity  limit  is  reduced  when  a null  s vec- 
tor solution  is  encountered  and  the  Feasible  Direction  problem 
reexamined  to  determine  if  the  point  is  indeed  an  optimum.  Since 
in  this  procedure 

Ej  - C2  otr/a°  (27) 

when  equation  (22)  is  satisfied  it  is  assumed  that  the  are 
sufficiently  small  to  allow  the  point  x£  to  be  considered  an 
optimum. 
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This  procedure  is  the  same  as  that  of  reference  [4]  except 
for  the  following  differences.  In  reference  [4]; 

1.  the  values  of  Vf(x^)  were  evaluated  at  each  point  calling 
for  the  computation  of  a penalty  parameter  X^ . Thus  Vf(x^) 
is  often  calculated  several  times  during  a local  explora- 
tion. Since,  however,  only  an  estimate  of  Vf  is  required 
here,  Vf(x^)  is  computed  here  only  once  during  local  explo- 
rations. A substantial  reduction  in  the  number  of  required 

% 

objective  function  evaluations  required  for  penalty  evalua- 
tions is  thereby  achieved. 

2.  all  constraints  were  evaluated  at  all  points.  Since,  how- 
ever, only  those  constraints  which  are  active  or  near  act- 
ive influence  movement  toward  the  optimum  at  a given  point, 
the  procedure  used  here  only  considers  such  constraints. 

All  constraints  are  evaluated  only  periodically  (after  each 
pattern  jump)  to  monitor  constraint  activity. 

3.  at  points  of  pattern  search  termination  in  the  unfeasible 
region  near  the  feasible-infeasible  boundary  the  pattern 
search  was  restarted  using  X^  - K,  the  large  penalty  para- 
meter, in  an  attempt  to  drive  the  search  to  the  feasible 
region.  In  the  strategy  used  here,  movement  in  a direction 
finding  problem  with  - 100,  is  utilized. 

4.  slightly  different  version  of  the  Rotating  Coordinate 
Pattern  Search  is  used.  The  version  employed  here  contains 
several  Improvements. 
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COMPARISON  STUDY 


The  work  o£  Eason  and  Fenton  [3]  forms  the  primary  basis 

for  the  comparison  between  the  algorithm  presented  here  and  other 

optimization  procedures.  A FORTRAN  IV  code  called  CAD0P3  based 

on  this  algorithm  was  developed  and  used  on  the  test  problems 

given  in  [3] . The  code  is  operational  with  IBM  FORTRAN  G and  H 

and  the  UNIVAC  TDOS  and  TSOS  systems.  This  code  required  the 

user  to  supply  only  the  objective  and  constraint  functions,  the 

% 

initial  values  of  the  variables,  and  regional  constraint  values, 
if  used.  The  convergence  criteria  and  initial  step  size  may 
also  be  specified. 

If  no  step  size  is  specified,  the  size  is  internally  gen- 
erated. In  this  program  the  initial  step  size  is  selected  such 
that  at  the  starting  point  a step  of  size  in  the  Vf  direction 
produces  a one  percent  change  in  the  value  of  the  objective  func- 
tion.  The  internally  generated  step  size  is  constrained  so  that 
> 0.005.  The  minimum  step  size  is  defined  as  ctm^n  ■ a^/1000 
unless  specified  otherwise  by  the  user. 

Nondimenslonal  constraint  equations  of  the  form 


W ■ <Bj 

- V/IV 

< 0 when 

0 

and 

\*  °1 

or 

> (28) 

w ■ <l3 

- v/iLji 

0 when  L j 

0 

and 

- °J 

are  used.  Otherwise  a dimensional  form  is  used. 

Here  represents  the  constraint  "behavior"  and  and 
Lj  the  upper  and  lower  limits  on  this  behavior,  respectively. 
The  user  supplies  the  expressions  for  f(x^)  and  the  expressions 
or  parameters  defining  B^ , U ^ , and  Lj . 
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A constraint  given  as  0 ^bCx^  <_  A would  be  written  as 

two  constraint  equations.  For  example,  one  could  be  of  the  form 

of  the  first  of  equations  (28)  where  B,  - b(x  ) and  U,  - A.  The 

in  1 

second  would  be  of  the  form  g„  - - b (x  ) since  L . - 0. 

2 n i 

The  program  uses  ej^  - 0.0001,  e2  - 0.10  and  e3  - 1 x 10~7. 
Multiple  starting  points  can  be  trfed  in  a single  computer  run 
to  provide  additional  optimality  confirmation.  Such  a technique 
is  recommended  even  though  the  algorithm  presented  here  seems 
reasonably  reliable,  since  no  nonlinear  method  can  be  considered 
absolutely  reliable. 

All  ten  problems  of  reference  [3]  were  run  with  the  CAD0P3 
code  from  the  starting  points  specified  in  reference  [3] . The 
internally  generated  step  sizes  were  used  for  all  problems  and  no 
algorithm  control  constants  were  changed  during  the  study.  Thus, 
there  was  no  tuning %of  the  code  to  individual  problems.  All  prob- 

/U 

lems  were  run  on  an  IBM  370  model  using  the  "G"  level  compil- 
er and  double  precision.  Thus,  the  runs  closely  simulate  those 
of  the  Eason  and  Fenton  study  which  likewise  used  an  IBM  machine 
with  the  same  level  compiler  and  precision  specification.  Thus 
the  comparisons  can  be  considered  accurate.  . - 

Table  1 briefly  describes  the  codes  studied  in  reference 
[3]  with  the  addition  of  the  DSDA  code,  the  earlier  CAD0P2  code, 
and  the  CAD0P3  code  studied  here.  Additional  details  on  these 
codes  can  be  found  in  references  [4]  and  [11] . Table  2 presents 
results  of  the  performance  of  these  codes  on  the  ten  test  prob- 
lems in  the  form  of  a success  matrix.  The  numerical  entries  indi- 
cate the  normalized  time  required  for  solution  [3],  the  symbol 


14 


CAD0P3 


TABLE  1 


CODE  DESCRIPTIONS 


Description 


ADRANS  Random  search  followed  by  pattern  moves 


class* 


Rosenbrock  search  •- 

Davidon-Fletcher-Powell  with  numerical 
derivatives 

Fletcher-Reeves  conjugate  gradient  method 
with  secant  approximation  derivatives 

Davidon-Fletcher-Powell  with  secant  approxi' 
mation  derivatives 

Hook  & Jeeves  pattern  search 

Steepest  descent  method 

Grid  and  star  network  search,  with  shrink- 
age 

Davidon-Fletcher-Powell  with  retained  step 
length  information 

Simplex  search 

Modified  pattern  search 

Modified  pattern  search,  dome  strategy 

Modified  pattern  search,  ridge  strategy 

Random  search  with  shrinkage 

Pattern  search  followed  by  random  search 

Modified  pattern  search 

Modified  simplex  search 

Modified  pattern  search  followed  by' 
Mugele's  search 


CAD0P2  Direct  search-feasible  direction  algorithm 


CLIMB 

DAVID 


DFMCG 


DFMFP 


FMIND 

GRADA4 

GRID1 


MEMGRD 


NMSERS 

PATSH 

PATRNO 

PATRN1 

RANDOM 

SEEK1 

SEEK3 

SIMP LX 

DSDA 


Refined  direct  search-feasible  direction 
algorithm 


DS  - direct  search,  DSF  ■ direct  search  employing 
SUMT  strategy  and  penalty  function,  G - gradient  proce- 
dure, GF  - gradient  procedure  using  SUMT  strategy  and 
penalty  function.  AR  - area  reduction  method. 


PERFORMANCE  OF  OPTIMIZATION  CODES* 


**Excluding  regional  constraints. 


"P"  indicates  progress  toward  a solution,  and  a blank,  indicates 
failure.  Definitions  of  "solution"  and  "progress"  vary  for  each 


problem  and  are  given  in  reference  [11]  as  well  as  a detailed 
description  of  each  problem. 

The  data  from  Table  2 may  be  applied  to  a number  of  rating 
schemes  for  comparing  the  codes  tested.  Table  3,  presents  relative 

rankings  of  the  codes  using  some  of  the  rating  criteria  of  ref- 
erence [3].  The-rating  equations  used  are  as  follows: 

The  number  of  problems  solved  by  code  "a"  is  nfl  and 

(29) 


N - n + n’  /2 
a a a 


where  n * la  the  number  of  problems  in  which  a P rating  was 

wL 


achieved.  The  efficiency  ratings  are  given  by 

10 


10 

Z 

- P"1 


b t /min(t  ) 
ap  ap  ap 


/n£ 


(30) 


where  bfi  ■ 1 if  code  "a"  solved  problem  "p"  and  zero  otherwise. 


t is  the  normalized  time  required  for  solution  (cpu  time  divided 


by  cpu  time  required  to  run  a timing  standardization  program) 


[3],  and  min(ta^)  is  the  shortest  time  required  by  any  of  the 


codes  to  solve  problem  "p."  The  other  efficiency  criterion  is 

10 


— 10 
Z 

_ P-1 


b t /mean(t  ) 
ap  ap  ap 


/n: 


(31) 


where  mean  (t  ) is  the  average  time  required  by  the  codes  to 
ap 


solve  problem  "p."  Here  the  overall  rating  number  for  general- 
ity (reliability)  and  efficiency  (speed)  is  given  by 


10 


T - Z t 

a p-1  ap 


(32) 
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where  tfi  Is  set  equal  to  twice  the  time  used  by  the  slowest  code 
solving  a problem  "p"  that  code  "a"  could  not  solve.  This  pen- 
alty time  Is  used  to  penalize  code  unreliability.  Only  codes  that 
solved  half  or  more  of  the  problems  are  rated  for  efficiency. 

The  tables  presented  here  are  essentially  similar  to  those 
given  In  reference  [3] , except  that  the  DSDA , CAD0P2  and  CAD0P3 
codes  are  Included.  Thus,  both  the  rating  schemes  and  presenta- 
tion of  results  are  essentially  those  of  [3]  . 

Based  on  Its  performance  on  the  ten  problems  of  reference 
[3]  the  new  algorithm  appears  to  be  a relatively  fast  and  reli- 
able optimization  procedure.  Codes  based  on  this  method  were 
the  only  ones  solving  all  problems  attempted.  CAD0P3  Is  compa- 
rable in  speed  to  the  faster  codes  (NMSERS,  and  PATRN1) . It  was 
significantly  faster  than  average  in  all  problems  and  was  the 
fastest  code  solving  problems  1 and  2.  Compared  to  the  rela- 
tively reliable  codes  (n  > 8) , CAD0P3  was  faster  than  DSDA  on 

a — 

seven  of  nine  problems  solved  by  DSDA,  faster  than  CAD0P2  on 
seven  of  ten,  faster  than  PATSH  on  seven  of  nine,  and  faster  than 
SEEK3  and  SIMPLX  in  all  problems  solved  by  these  schemes.  As 
may  be  seen  that  DSDA  and  CAD0P2  are  the  only  reasonably  reli- 
able codes  comparable  in  speed  to  CAD0P2.  CAD0P3  appears  signif- 
icantly faster  than  PATSCH,  SEEK3,  and  SIMPLX.  Only  the  straight 
(PATRN1)  and  simplex  (NMSERS)  codes  were  generally  faster.  Of 
this  group,  only  NMSERS  appears  to  be  sufficiently  reliable  to 
merit  serious  consideration.  The  differences  in  speed  between 


NMSERS  and  CAD0P3  is,  however,  overshadowed  by  the  superior  reli- 
ability of  the  latter.  Thus,  in  the  overall  generality  and  effi- 
ciency rating,  CAD0P3  stands  alone.  Viewed  on  the  basis  of  these 
comparisons,  CAD0P3  appears  to  be  a superior  nonlinear  MP  code. 

It  should  be  noted  that  problems  4,  5,  8 and  10  are  uncon- 
strained at  the  optimum  points.  Thhs,  even  though  problem  10  has 
a regional  constraint  which  is  active  at  the  start  of  the  search, 
these  problems  can  be  considered  essentially  unconstrained. 

v _ 

Those  codes  (PATRN1  and  (NMSER)  which  are  faster  overall  than 
CAD0P3  obtain  their  superior  efficiency  ratings  primarily  from 
faster  solutions  of  such  problems.  The  DSFD  algorithm  is  de- 
signed, however,  for  constrained  problems.  Thus,  although  its 
overall  performance  on  constrained  and  unconstrained  problems, 
or  its  performance  on  unconstrained  problems,  is  of  interest,  a 
direct  evaluation  of^  its  performance  on  constrained  problems  is 
also  appropriate.  Thus  a comparison  on  those  problems  with  be- 
havior constraints  (problems  1,  2,  3,  6,  7,  and  9)  is  given  in 
Table  4. 

In  Table  4 only  those  codes  that  solved  three  or  more  of 
the  constrained  problems  are  evaluated  for  efficiency.  In  this 
comparison  CAD0P3  is  comparable  to  NMSERS  in  speed  with  these 
codes  being  substantially  more  efficient  than  any  other  code. 
CAD0P3  is  ranked  much  higher  in  the  overall  rating  since  NMSERS 
failed  on  two  of  the  six  problems.  CAD0P3  is  faster  than  CAD0P2 
and  DSDA  on  all  but  one  problem.  This  problem  has  only  one  con- 
straint, however,  and  thus  the  reduction  associated  with  checking 
only  nearby  rather  than  all  constraints  is  lost. 
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TABLE  4 


RELATIVE  RANKING  OF  OPTIMIZATION  CODES  ON  PROBLEMS 
WITH  BEHAVIOR  CONSTRAINTS 


Generalty 


Efficiency 


Generalty 

and 


n 

N 

f 

f 

T 

a 

a 

a 

a 

a 

& 

CADOP3 

CADOP3 

2.2 

NMSERS 

0.14 

NMSERS 

1.2 

CADOP3 

CADOP2 

G CAD0P2 

2.6 

CADOP3 

0.16 

CADOP3 

2.0 

CADOP2 

DSDA 

DSDA 

6.3 

CADOP2 

0.30 

CADOP2 

8.1 

DSDA 

PATSH 

. PATSH 

J 

ADRANS 

5 ADRANS 

12. 

DSDA 

0.44 

DSDA 

8.8 

NMSERS 

SEEK3 

5 SEEKS 

22 

PATSH 

0.93 

PATSH 

11 

13 

PATSH 

4 

NMSERS 

4.5  NMSERS 

SEEK3 

1.1 

NEMGRD 

DAVID 

SEEK3 

SIMPLX 

_ _ SIMPLX 
3 MEMGRD 

25 

NEMGRD 

NEMGRD 

MEMGRD 

27 

DAVID 

1.2 

SEEK3 

14 

ADRAMS 

3 

SIMPLX 

DAVID 

DAVID 

51 

SIMPLX 

1.9 

SIMPLX 

18 

DAVID 

- 0 FMIND 
SEEK1 

65 

ADRANS 

2.9 

ADRANS 

FMIND 

RANDOM 

PATRN1 

SEEK1 


PATRN1 

2.5  RANDOM 

_ _ PATRNO 
* GRID1 

1.5  GRAD4 


PATRNO 

GRIDI 

GRAD4 

DFMCG 

DFMFD 

CLIMB 


1.0  DFMFD 
0.5  DFMCG 
0.0  CLIMB 
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CAD0P3  is  faster  than  NMSERS  on  three  of  four  problems  although 
both  codes  are  essentially  similar  in  speed.  It  is  faster  on 
all  problems  than  the  remainder  of  the  codes  rated  for  efficiency. 
CAD0P3  is  first  in  the  generality  and  efficiency  rating  since  it 
is  faster  than  CAD0P2  and  since  all  other  codes  were  penalized 
in  this  rating  for  failure  on  one  or  more  problems.  It  may  be 
seen  therefore  that  the  DSFD  improvements  cited  here  did  indeed 
increase  performance. 

Thus  the  improved  DSFD  procedure  seems,  on  the  basis  of 
this  comparison,  clearly  superior  to  the  other  procedures  studied 
on  constrained  problems. 

Table  5 lists  the  number  of  function  evaluations  required 
by  CAD0P3  to  reach  solution.  The  comparison  data  of  reference 
[3]  does  not  provide  these  values  for  the  other  codes.  These 
values  are  provided  here  since  it  is  felt  that  for  most  engineer- 
ing problems  the  number  of  function  evaluations  is  a better  cri- 
terion of  performance  than  nondimenslonal  time. 

The  time  required  for  solution  can  be  viewed  as  the  sum 
of  the  function  evaluation  and  algorithm  execution  time.  For 
very  simple  functions  most  of  the  required  time  may  be  associated 
with  algorithm  execution  particularly  where  the  algorithms  are 
complex.  For  such  problems,  however,  time  is  usually  of  little 
Importance  since  any  reasonably  effective  procedure  can  generate 
a relatively  low  cost  solution.  On  the  other  hand,  efficiency 
becomes  important  where  function  evaluation  Is  computationally 
demanding.  Here  algorithm  execution  time  may  be  quite  small 
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TABLE  5 


NUMBER  OF  FUNCTION  EVALUATIONS  FOR  SOLUTION 


Constraint 

Objective  Function  Function 

Problem  No. Evaluations Evaluations 


1 

2088 

4860 

2 

266 

190 

3 

622 

993 

4 

4910 

— 

5 

347 

— 

6 

378 

322 

7 

556 

737 

8 

217 

— 

9 

488 

447 

10 

849 
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compared  to  function  evaluation  time.  In  such  cases  the  number 
of  function  evaluations  is  a reasonably  accurate  criterion  for 
performance . 


Preliminary  experiments  indicate  that  although  nondimen- 
sional  time  is  a valid  comparison  criterion  where  codes  are  run 
using  the  similar  compilation  software  it  is  not  valid  otherwise. 
It  is  felt  therefore  that  the  number  of  function  evaluations  is 
a superior  comparison  criterion  since  it  is  a reliable  indicator 
of  performance  on  problems  where  efficiency  is  important.  Thus 
this  information  is  provided  here  so  as  to  allow  future  compar- 
isons . 

The  principal  apparent  disadvantage  of  CAD0P3  is  its  rela- 
tive complexity  compared  to  some  of  the  other  reasonably  reliable 
methods.  It  contains  713  FORTRAN  statements,  while,  for  example, 
DSDA  contains  372  and  PATSH  only  about  75.  Thus  for  simple  prob- 
lems of  relatively  low  dimensionality,  one  of  the  simpler  codes 
may  be  preferable. 

It  may  be  seen  that  the  gradient-based  procedures,  includ- 
ing those  employing  the  SUMT  strategy  and  penalty  function  [12], 
so  widely  used  in  engineering  problems,  performed  rather  poorly. 
The  best  gradient-based  code  DAVID  solved  only  half  the  problems. 
The  direct  search  procedures  in  general  provided  much  better  per- 
formance although  only  four,  CAD0P3,  CAD0P2,  DSDA,  and  PATSH, 
solved  or  approached  the  solution  to  all  ten  problems  with  only 
CAD0P3  and  CAD0P2  solving  all  problems. 
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It  should  be  noted,  however,  that  this  comparison  employed 
rather  small  test  problems  (two  to  five  variables  with  zero  to 
ten  behavior  constraints) . It  is  not  clear  that  the  direct  search 
procedures  would  also  possess  reliability  and  speed  superior  to 
the  gradient  methods  on  large  problems.  Still,  the  superiority 
in  this  study  of  the  better  direct  search  procedures  is  dramatic 
while  the  performance  of  the  gradient  schemes  is  rather  dismal. 
Furthermore,  although  the  better  direct  search  procedures  may  not 
be  as  reliable  on  larger  problems,  one  would  certainly  not  expect 
the  reliability  of  the  gradient  procedures  to  Improve  on  such 
problems . 
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CONCLUSION 


Although  it  Is  recognized  that  such  favorable  findings  on 
the  CAD0P3  code  presented  by  Its  developer  can  be  viewed  with 
some  skepticism.  It  should  be  noted  that  the  performance  cited 
Is  based  primarily  on  rating  schemes  and  problems  selected  by 
Eason  and  Fenton  and  not  by  the  author.  Furthermore,  this  per- 
formance was  achieved  with  a conscious  effort  to  avoid  any  tuning 
of  the  code  to  Individual  problems.  Thus,  these  results  Indicate 
that  the  DSFD  algorithm  appears  to  be  a superior  nonlinear  mathe- 
matical programming  procedure  at  least  on  relatively  small  problems 
on  which  Is  was  tested. 

It  should  be  noted,  however,  that  this  algorithm  combines 
the  strength  of  the  direct  search  procedures  (speed)  and  the 
Zoutendljk  method  (ability  to  confirm  optimality  and  identify  a 
feasible  direction)  In  such  a way  that  the  weakness  of  the  direct 
search  In  treating  problems  of  large  dimensionality  is  minimized. 
Thus,  one  would  expect  this  new  procedure  to  also  be  quite  effec- 
tive on  larger  problems.  This  potential  h.as,  however,  yet  to  be 
demonstrated . 
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APPENDIX  A 


USER  INSTRUCTIONS  FOR  CADOP3 

1 . Instruction 

This  code  treates  inequality  constrained  or  unconstrained 
linear  or  nonlinear  minimization  problems  by  means  of  a direct 
search  mathematical  programming  method.  The  method  starts  from 
a user  specified  starting  or  initial  point  and  generates  a se- 
quence of  better  points  until,  hopefully,  an  optimal,  or  near 
optimal,  point  is  reached.  The  code  is  intended  primarily  for 
constrained  nonlinear  problems.  Linear  problems  are  much  more 
effectively  handled  by  one  of  many  linear  programming  methods. 
Unconstrained  problems  are  better  treated  by  UNCDP,  a simpli- 
fied version  of  CAD0P3,  or  one  of  several  other  unconstrained 
nonlinear  problem  codes.  CAD0P3  will,  however,  treat  both 
linear  and  unconstrained  problems. 

It  should  be  noted  that  nonlinear  mathematical  programming 
methods  generally  do  not  guarantee  an  optimal  solution  because 
of  the  possibility  of  multiple  local  optima,  algorithm  failure, 
or  numerical  difficulties.  Alghough  CAD0P3  is  among  the  most 
reliable  of  the  nonlinear  mathematical  programming  codes  it 
cannot  be  considered  absolutely  reliable.  Thus,  it  is  desir- 
able to  use  several  optimization  runs  with  widely  separated 
starting  points  to  confirm  the  convergence  has  been  achieved 
or  to  determine  if  local  optima  are  present. 
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2 . Problem  Formulation 


This  code  treats  the  problem: 

Find  those  values  x^^  of  the  set  of  "variables"  x^  i ■ 1,2 IP 

and  the  constant  "parameters"  p^,  k *»  1,2 KP,  that  result  in 

the  minimum  value  of  the  objective  function  f (P^*3^)*  that  is 

f (pk.**iP  “ mln  Cl) 

while  also  satisfying  the  "behavior"  constraints 

gj  - 0 J - 1>2» JP  (2) 

and  the  "regional"  constraints 

xi  1 xi  5 xi  (3) 

Two  types  of  behavior  constraints  are  recognized,  one 
type  are  "linear"  constraints  of  the  form 

IP 

(,1>  ' i-i  ‘ii*1  <4) 

Where  f (x^)  is  also  of  this  form  one  has  a linear  programming 
problem.  If  any  of  Eqs.  (1  and  2)  is  not  of  this  form  the 
problem  is  nonlinear.  Constraints  not  of  the  form  of  Eq.  (4) 
are  nonlinear  constraints. 

The  program  in  its  present  form  treats  problems  with 
- 1 Ip  1 10»  0 1 JP  1 1°  and  < KP  < 100 


- 30  - 


... 


■ T' 


3 . Specification  of  Oblectlve  and  Constraint  Functions 


The  subprogram  FUNCTION  OBJ(J)  is  essentially  a dummy 
subprogram  created  to  accept  the  user s FORTRAN  IV  program 
statements  defining  the  objective  and  behavior  constraint 
functions.  The  normal  form  of  the  behavior  constraints  is 


B1  (pk,xl>  - Uj  (pk»xiP 


(4) 


'j  - i 

where  B^  can  be  thought  of  as  the  behavior  and  the  upper 
limit  of  behavior. 

The  FUNCTION  OBJ(J)  is  modified  so  as  to  define  the 
objective  and  constraint  functions  by  inserting  the  following 
statements  immediately  after  the  COMMON  statement 

GO  TO  (1,2, jp),  J (omit  if  no  constraints  are 

used) 

OBJ  * expression  defining  f using  arrays 

P(k),  X(i) 

RETURN 

j B = expression  defining  B 


U * expression  defining  U 
GO  TO  101 


A constraint  statement  group  is  used  for  every  constraint. 

The  last  constraint  group  need  not  use  the  final  GO  TO  state' 
ment . 
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For  the  problem 
minimize 


i 

> 

\ 


l 

f 

I 


Example 


2xl  X2  " P1X3X4 

such  that 

2x33  - X1  - P2 

and 

2 

“6  <.  X2  - x^  <_  0 

one  could  define  three  constraints  and  insert  the  statements 
GO  TO  (1,2,3) ,J 

OBJ  - 2.  * X (1)  * X (2)  - P (1)  * X (3)  * X (4)  **2 
RETURN 

1 B - 2.*  X (3)  **  3 - X (1) 

U - P (2) 

GO  TO  101 

2 B - X(2)  **  2 - X (4) 

U - 0. 

GO  TO  101 

3 U - X (2)  **  2 - X (4) 

B - - 6. 
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Alternately,  noting  that  both  upper  and  lower  limits 
on  the  second  constraint  equation  cannot  be  active  simultane- 
ously, one  could  use  two  constraint  equations  and  insert  the 
program  segment: 


GO  TO  (1 , 2)  , J 

OBJ  - 2.  * X (1)  * (2)  - P (1)  * X (3)  * X (4)  **  2 
RETURN 

l 

1 B - 2.  * X (3)  **  3 - X (1) 

U - P 
GO  TO  101 

2 TB  - X (2)  **  2 - X (4) 

IF (TB . LT , -3  . ) GO  TO  3 
B - TB 

U - 0. 

GO  TO  101 

3 B - -6. 

U - TB 

The  constraints  values  at  the  optimum  are  printed  in 
the  form  of  equation  (2)  where 


(Bi  - 

(Bj  - 0j) 


Uj  4 0 and  B ^ 4 0 
Uj  - 0 or  B^  - 0 


thus  a negative  value  of  g^  indicates  the  constraint  is  sat- 
isfied . 


■rmt 

j 


33 


• '1 

| i 


A ■ 


r 


4 . Data  Cards 

The  program  allows  multiple  optimization  runs  of  a par- 
ticular problem  within  a single  computer  run  using  different 
parameters,  for  parametric  studies,  different  starting  points, 
for  optimality  confirmation  or  search  for  local  optima,  and 
different  sets  of  regional  constraint  limits. 

The  user  must,  by  use  of  the  data  card  set,  specify  the 
number  of  parameters,  variables  and  behavior  constraints  used, 
identify  the  linear  constraints,  specify  for  the  first  run  if 
regional  constraints  are  to  be  used,  and  for  subsequent  runs 
whether  the  parameters.  Initial  point  or  regional  limits  are 
to  be  changed.  If  parameters  are  used  they  must  be  entered. 

The  initial  point  must  be  entered.  The  values  for  the  regional 
constraints  must  be  entered  if  such  limits  are  imposed. 

The  first  data  set  is  entered  on  the  "Problem  Constants" 
card  which  specifies.  In  order,  the  number  of  parameters  (KP), 
variables  (IP)  and  behavior  constraints  (JP)  , used  (see  line  9 
of  program  listing).  The  data  is  entered  in  3110  format  in  the 
first  3 fields  in  order  and  must  be  right  justified. 

If  constraints  are  used  (JP  > 0)  a second  data  set 
is  entered  on  the  "Constraint  Identification"  card  which  iden- 
tifies the  linear  behavior  constraints.  If  a constraint  is 
linear  a digit  is  entered  in  the  column  with  the  same  number 
as  the  constraint.  If  no  linear  behavior  constraints  are  used 
a blank  card  is  employed  (see  line  10) . 


The  third  data  set  is  entered  on  the  "Data  Control"  card 
and  is  used  to  specify  whether:  1)  new  parameters  are  to  be 

used  (ICNTR2),  2)  a new  initial  point  is  to  be  used  (ICNTR3) 
and,  3)  regional  limits  or  new  regional  limits  are  to  be  used 
(ICNTR4)  (see  lines  26-43)  . The  entries  are  made  in  3110  format 
in  the  order  given  above.  For  the  first  run  no  entries  in  the 
first  two  fields  are  needed.  If  regional  limits  are  used  an 
entry,  any  nonzero  digit,  is  made  in  the  third  field  (col  21-30). 

If  and  only  if  the  number  of  parameters  specified  is 
greater  than  zero  (KP  > 0)  a fourth  data  set  is  used  to  enter 
the  problem  parameters  p^.  The  entries  are  made  5 to  a card 
in  F15.5  format  in  order  i ■ 1,2 KP , (see  lines  21  and  22). 

The  fifth  data  set  specifies  the  starting  point  and 
entries  are  made,  in  order,  8 to  a card  in  F10.5  format  (see 
lines  23  and  29) . 

If  and  only  if  an  entry  is  made  in  the  third  field  of 
the  Data  Control  card  (ICNTR4  ■ 0)  a sixth  data  set  defining 
the  lower  limits  is  entered,  8 to  a card,  in  F10.5  format  fol- 
lowed by  an  upper  limit  set  (starting  with  a new  card)  with 
similar  format.  Both  complete  limit  sets  must  be  entered  (see 
lines  30-34 ) . 

For  every  additional  run  an  additional  Data  Control  card 
is  added  followed  by  parmeter  and/or  initial  point  and/or  re- 
gional limit  sets,  where  and  only  where  the  need  for  a new  set 
is  indicated  by  an  appropriate  entry  (ICNTR2  - 1 or  (ICNTR3  * 1 
or  ICNTRL4  ■ 1)  on  the  Data  Control  card  (see  lines  26-34). 

A blank  Data  Control  card  is  used  to  terminate  the  run(s) 
(see  lines  26  and  27). 
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5 . Conversion  To  Single  Precision 


The  program  as  supplied  is  a double  precision  form. 
Experience  has  shown,  however,  that  the  great  majority  of 
problems  can  be  treated  adequately  in  single  precision  and 
thus  it  is  recommended  that  single  precision  be  tried  first. 
To  convert  to  single  precision  remove  the  IMPLICIT  statement 
at  the  beginning  of  MAIN  and  all  subprograms  and  the  function 
conversion  statements  such  as  SQRT(B)  * DSQRT(B)  etc.,  near 
the  beginning  of  most  subprograms.  Furthermore,  the  REAL 
FUNCTION  0BJ*8(J)  should  be  put  in  single  precision  form. 
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6 . Change  of  Problem  Size 


To  change  the  maximum  number  of  variables  (IP)  or  the 
maximum  number  of  constraints  CJP)  the-  program  can  treat, 
change  the  arrays  in  MAIN  and  all  subprograms  as  follows: 


1.  change  all  D,  DZ,  X,  DL , and  DU  arrays  in  common  to 
D (IP ) , DZ (IP)  etc., 

2.  change  the  G array  (T  in  0PT2  and  LINPRO)  as  well 
as  the  CK,  1A,  IC  array  in  COMMON  and  the  BL  array 
in  COMMON/OPT/  where  they  occur  to  G(JP),  CK(JP) 

ETC.  , 

3.  change  the  V array  (A  array  in  0PT2  and  LINPRO)  to 
V (M ,N)  in  COMMON  where 

M = JP  + 4IP  + 5 
N - JP  + 5 IP  + 6 

4.  change  all  DUMMY  arrays  such  as  DUM,  IDUM,  etc., 
to  equalize  common  sizes. 

5.  In  the  0PT2  SUBROUTINE  DIMENSION  statement  change 
DT,  SS,  IB,  IC,  to  DT(IP)  etc.,  change  GG  to  GG(JP) , 
TN  to  TN ( JP+1) , ID  to  ID ( JP) , DK  to  DK(JP+1,  IP), 

DBJ  to  DBJ  (JP+1,  IP)  and  SSS  to  SSS(M)  where  M 

is  as  above. 


6.  In  the  FIND  SUBROUTINE  DIMENSION  statement  change 
DT (10)  to  DT(IP) . 


7.  In  the  PATTRN  SUBROUTINE  DIMENSION  statement  change 
all  10's  to  IP. 


MAIN 


IMPLICIT  REAL*8(A-H,0-Z) 

CCMMON  D110),P(1G0),0Z(1QI,X(10),C  (10  ),A,AM IN,DL ( 1C), DU (10)  • 
1VI55,66),CK(10>,Q,IP,JP,IA(10» ,KP  ,KW,KL, IK, I0UMI21I tJX 
COMMON/QP  T/BL ( 10 1 ,TM 
1 FORMAT ( 5F 15 .5  I 

2 FORMAT (4110) 

3 FORMAT (2F10. 51 

101  FORMAT (8F1Q.5  1 

READ2, KP, IP  yJP 

1FUP.NE.0IREA0  44,  (CK  (K ) ,K  = 1 , JP  I 
1FUP.EO.OIGO  TO  45 
00  46  KS1 , JP 

IF(CKtK).NE.0.)CK(Kl*0.01 
1F(CK(K!.EQ.0.)CK(KI*1  . 

1F(CK(K).LT..5)CK(K}*0. 

46  CONTINUE 
45  OE  21  1*1, IP 
DL(I )=-l.E*48 
44  FORMAT (80F1.0  I 
21  DU(I)*1.E*48 

READ2»1CNTR2,ICNTR3,ICNTR4 
1 F (KP.EQ.O  1 GO  TO  12 
READ1, IP(11,1*1,KP1 
12  READ101,(D(I1,I*1,1P> 

18  FORMAT ( 1H0 , • STARTING  DESIGN  VARIABLE  VALUES*  I 
GO  TO  23 

24  READ(5,2,END*33)ICNTR2«ICNTR3, ICNTR4 

IFC1CNTR2+ICNTR3+ICNTR4.EC.0IST0P 
IF(ICNTR2.NE.0)READ1,  (P  ( I I , I *1  ,KP  ) 
IFIICNTR3.NE.0)READ1Q1,(D( I)  ,1*1,  IP) 

23  1FOCNTR4.EQ.OIGQ  TO  22 
READ  101,(DL(I),1*1,1P) 

17  FORMAT ( 1H0,  • LOWER  LIMITS  OF  DESIGN  VARIABLES*) 

READ  101»(DU(I),I*1,IP) 

19  FORMAT ( 1H0 , * UPPER  LIMITS  OF  DESIGN  VARIABLES*) 

22  IF(KP.LE.O)  GO  TQ  104 

PRINT  6 

PRINT7,(K,P(K),K*1,KP) 

104  PRINT  18 

PRINT 101,(0(1), 1*1, IP) 

PRINT  17 

PRINT 101, (DL(II,IX1, IP) 

PRINT  19 

PRINT 101, (DU( I ) , I - 1 , 1 P ) 

CALL  PATTRN 
IK*0 

CALL  FIND(TO) 

OB JD*OB J (0  I 
PR INT4 ,KL 

4 FORMAT ( *0N0.  OF  REDESIGN  CYCLES**, 15) 

PR1NT5 ,OBJD 

5 FORMAT  ( 'OOPTZMUM  VALUE  * * ,F 13  .8  ) 

6 FORMAT!  *ODESIGN  PARAMETERS*/) 

7 FORMAT l • P * ,1 2 » •* • ,F14  »4) 

PRINT8 

8 FORMAT ( 'ODES 1GN  VARIABLE  VALUES*/) 
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MAIN 


PRINT9,  (K,D(K1 ,K>1,IP) 

9 FORMAT t • OSU.'^SFlS.e) 

IFUP.EQ.OI  GO  TO  24 
PRINTIO 

10  FORMAT ( 'ONEARNESS  TO  CONSTRAINTS  •/ I 
PRINT  11, {K,G(K1 ,K*1,JP ) 

11  FORMAT  ( • GSU.^SFIS.S) 

100  GO  TO  24 

33  STOP 

END 


>•**«*;* 


PATTRN 


SUBROUTINE  PATTRN 

IMPLICIT  REAL*8(A-H,C-Z1 

CCMHON  DtlQl,P(1001,DZ(10!  ,X (10  I ,6 (10  I , A, AM  IN ,OL 1 10 ) ,DU f 10  1 
lVf55,661,CK|10),0,IP,JP,IA(10},KP.KN,KL,lK,IDUMI2Ii,JX 
COMMON/OPT/BLdO ) ,TN 

C1MENSI0N  AO (10, 10 >,E  1(10,10 1 ,ET (10,10  I ,AZ( 10  I 
TCIF-10.E-U0 
IPP*IP 
AL*A 
111*0 

STE5T*0. 

KA*0 

KU*0 

KT*0 

KL*0 

KF*0 

KT*0 

KWN*0 

CG  112  K*1 , JP 
112  BL (Kl *-0 .1 

DC  111  K*1 , I PP 

IFIDtK  ).LT.0L(KH0(K1*0L(K  ) 

IF  (0(K  ) .GT.DU  (KUO  (K  I*OU(K  I 
111  X (K  1*0  (K I 
CALL  SIZE 
25  KK=  1 

DC  1 K*1 , 1 PP 
00  22  I*1,1PP 
EI(I»K)*0. 

ETU.K  1=0. 

22  ADCI,KJ*0. 

ETIK.K 1=1. 

E 1 (K,K 1*1. 

1 CONTINUE 
5T*SB 

1F(KA.EQ.11G0  TO  9 
KA*1 
190  1K*0 

JX*0 

CALL  FIN0(SB1 
9 ST*SB 

JX*0 

GO  TO  100 
19  1K*0 

JX*0 

CALL  FIND (ST  I 
100  KL*KL*1 
IK*1 
KF*KF-»i 
DO  2 J*1,IPP 
DC  17  K*1 , 1 PP 
AZ(K>*  E I ( J,K I *A 
1F(KK.NE.01AZ(K1«ETU,KI*A 
17  D(K1*0(K)4AZ(K) 

DC  35  K*1,1PP 

1F(0(K  I .GT.OUtMlGO  TO  4 
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it: 
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1F(D(K).LT.DL(K))GG  IQ  4 

35  CONTINUE 
CALL  FIND(SA) 

IF( (ST-SA1/DABS(SBI.GT  .2.E-16)  CO  TO  5 
6 DO  18  K*1 , l PP 
18  DCKI»DIK>-AZ(K»-AZ(KI 
OC  36  K*1 * l PP 
IF(D(K1  .LT.DLIKMGC  TO  3 
1F(0(K).GT.DU(K!)GO  TO  3 

36  CONTINUE 
CALL  FINO(SA) 

1 F ( (ST -5 A ) /OABS (SB  1 «LT  .2 .E-lfc I GO  TO  3 
5 ST*SA 
GO  TO  2 

3 DC  21  K*1 * I PP 
21  D(Kl«D(KI*AZ(K) 

2 CONTINUE 

1 FC ISB-ST1/DABS(SB».LT.2.E-16)  GO  TO  7 
26  SUM»0. 

DO  8 K *1  * 1PP 
DB«2.*D(K1-X(K) 

AD(1,K)*D(K)-X(K1 

1 F (DABS (AD ( 1, K > I. LT. *1.0-321  AD(1,K1*0 
SLM*SUM+ADll,K J*AD(1,K) 

X CK)*D ( K) 

S B =5T 

IFIDB.LT.DLfK  ) )DB=DL (K  ) 

IF(DB.GT.0U (K  I 1DB-CU  CK ) 

8 D(K)-DB 

SUM  =D$QRT (SUM  I 
IFISUK.LE.OIGO  TQ  23 
IFCSUH-DABSI  .5«AM  7,7,950 
950  CONTINUE 
KK*0 
SUH=0 ■ 

DC  11  J*1,1PP 
ET(1,J)*EI (1,J) 

E1(1,J)*AD(1,J) 

11  SUM*SUH+AD(1, JI*«2 
SUM*DSQRT (SUM  I 
DO  53  J*1,IPP 
53  E I(1,J)*EK1,J)/SUM 
DO  13  I *2 , I PP 
K*  1-1 

00  16  J*I,IPP 
16  AC ( I, J 1*AD (1, J I 
SUMA-O. 

DO  15  J-l.IPP 

IF(0ABS(E!(K,J)).LT.«l.E-32)  EI(K,JI«0 

15  SUMA*SUMA+E1 (K  , Jl*ADi  I , J) 

SUH»0. 

DO  16  J«1,1PP 

ETC  I , J 1«E I ( I, J 1 

EI(1,J)*AD(I,J]-EI(K,J)«SUMA 

IF (OAB  S(E1(I,JH.LT.+1 .E-32 ) EI(1,J)*0 

16  SUH«SUM«E1 (I,J)»*2 
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SUM*DSCRT I SUM  1 
IFiSUH )27,27,28 

27  DC  29  J*1,IPP 

29  E 1 C I » J )*ET  1 1 1 J ) 

GC  TD  19 

28  DC  13  J«1,IPP 

13  EI(1»JJ*EICI»J1/SUH 
GC  TO  19 
7 DC  20  K*i # I PP 
20  D(K)«X(K) 

1FCKKI25, 24,23 
24  KK*-1 

I F (KF .GT.200 1G0  TO  23 
€0  TO  9 
23  CONTINUE 

TN*5B 

IFIJP.EQ.OJGO  TO  32 
IK*0 

CALL  FIND ( SB  1 
DO  30  l«l,JP 

1 F C G C I ) .GT.-.01*BL  Cl  I) GO  TO  31 

30  CONTINUE 

IFCKKW.LE.3J  KH«0 
Kt«0 

GC  TO  32 
31  Kb*KW«l 

KT«KV*1 
KWb*KWW+l 
IK*0 

IFCKH.NE.O)  CALL  FINOCSBJ 
32  IF (KL .GT .KT*4  ) GO  TO  33 

A*A/2. 

00  34  K*1 » JP 
34  BLIKl *BL (K  J/2. 

33  CALL  0PT2C5B»ST»JK»KYJ 
KT«KL 

IF ( JK.EQ.2 I RETURN 
311  KF«0 

CO  300  1*1, JP 

IFCGCl ) .GT .-.0 1*BL  C I ) ) GO  TO  301 

300  CONTINUE 
IFCKWW.LE.3J  KW*0 

Kt*0 

301  IFCJK.NE.O.J  GO  TO  46 
IFCA.EQ.ALI  GO  TO  26 

PRINT  114,A,SB,ST,CDtl!,I*l,IP) 

114  FORMAT ( * 0PT2BA  SEA*  •»F11.9, *SB*s,E16.8t “ST *,»E16.8,,X*,,10F10.6J 

AL*A 

DIF*DABSt  CSB-STEST1/SB J 
IFCDIF.GT ..000011  GO  TO  45 
IP»1H*1 

IFCIM.LT. 2)  GO  TO  150 
IFCDIF.GT. TDIF J GO  TO  45 
RETURN 
45  1M*0 

149  STEST*SB 
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148 


113 

46 

47 


150 


115 


TO  I F*D  IF 
GO  TO  26 
IM*0 
l*.5*A 

OO  113  K*1,JP 
BL (K 1 «BL IK  1/2 • 

OC  47  K*1 * I PP 
0 IK  ) =X l K 1 

IF(A.LT.AMIN)  RETURN 
GO  TO  190 
A *A/2 

CO  115  K*1,JP 
8L(Kl =BL (K 1/2. 
IFtA.LE.AHINl  RETURN 
GO  TO  149 
EftD 


|i 

I* 

i -h 
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FIND 


SUBROUTINE  FINDITOl 

IMPLICIT  REAL*8  U-H.O-Z) 

CCMMON  D(10»,P(lQ01,OZ1101,XtlQI,G(10),A,AHIN,DL(10),DUC10l, 

IV (55.66) ,CK(10),Q,IP,JP,IA(10>,KP,KW,KL,IK,!L.IDUM ( 10 1 , IC ( 10 ) , JX 
COMMON/QPT/BL (10  I.TN 
C 1 HEN  5 ION  DT(IO) 

1*0 

TC-OBJ(O) 

TN*TO 

1FUP.EQ.0)  RETURN 
OC  25  K*1,IP 
25  D1(K1*D(KI 

IF(IK.EQ.O)  CO  TO  51 
IF(IL.EQ.O)  RETURN 
CO  52  J*1,1L 
I8«IC  U) 

G( IB) *0B  J (IB  I 
IF (G ( IB  ) I 52.52,5 

5 1*1*1 
IA( 1 1 *1B 

52  CONTINUE 

GO  TO  56 
51  11*0 

OC  3 J*1,JP 
G ( J) =0B J( J ) 

IF(G(J).LT.2.«BL(J))G0  TO  3 

6 IL*IL*1 
1C(IL )*J 

IF(GCJ))  3.3.4 
4 1*1*1 
I A ( I )*J 
3 CCNT1NUE 
56  I F ( I )4 1 .41.100 

100  fcTEHP*0. 

00  21  K *1 , 1 
IB* IA  (K 1 

IF(KU.NE.O)  GO  TO  12 
1FIG(IB).L1.-BI(1B>)G0  70  21 
72  E*10000.*G(  IB) 

IF(E-GT.HTEMP)  HTEMP*E 
21  CCNTINUE 

1F(HTEMP.G1..001)  GO  TO  19 
IF(JX.NE.O)  GO  TO  42 
SUM*0.0 
JX*1 

00  1 K *1 , IP 
DT(K)-D(K) 


FIND 


TA-0BJ101 
DC  18  K*1 , 1 
1B*IA(K ) 

TT*08J ( |B  I 

If  (TT-GCIBM  37,54,37 
37  e>DABS((2.*(T0-TA)y(lT-G(IBm*GUBn 
1 F (G l IB 1*BL( IB  1) 55,55,54 
54  E *10000 .«G { IB  I 
55  IF(E.GT.HIEHP)  HTEMP«E 

18  CCN11NUE 
19  CONTINUE 

TO=TQ+HTENP 
OC  26  K*1 ,IP 
26  DIKJ*DT(K1 
41  RETURN 
END 
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SUBROUTINE  0PT2 (SB ,ST » JK ,K Y l 
IMPLICIT  REAL«8 (A-H.O-2) 

INTEGER  Z ,E,G,BB,W,B,H 

CCHMON  D ( 10 ) ,F (1001  ,DZ (10 ) »X  1 10  ) »T  (10  ) ,V,AM1M,DL(  IC  I ,DU|  101  , 
1*155,66  I, CM  10  I ,Q,lP,JP,IA(lOl,KP,KW,KL,lK,IL,M,N,B,L,BB,E,G,Z  ,H  , 
2M  , 10UM  ( 10  1 • J5 

COHMON/QPT/BL (10  I ,TP 

DIMENSION  DT(10),SS(101,IB(10I,IC(10I,SSS(55),GG(1C),TN(11I,0K(101 
1,201, ID(1CI,DBJ(11, 10} 

VT*V/1 0 • 

TQ*TP 

J6*0 

JIN*0 

JXN*0 

CEL*. 00001 
IF(V.LT.DEL)  DEL*V 
IPI*JP*A«IP+A 
CO  116  K*1 , JP 
CG (K I *T (K  1 

116  ID (K } *0 

117  CO  103  1*1,  IPi 

103  SSS(I1*0. 

JK*1 

N*2*IP*1 .0 
N*IP*1 
Z*-l 
MM*JP«1 
M*HM+2«N*2«IP 
L*MM-*N*  IP 
E*0 
C*0 
NN*6 

90  0*99999 

190  B*M*N*G*1 

N*M 
320  H*1 

BB*B+1 

J*0 

DC  A 1*1, IP 
IB (1 1*0 

icm*o 

1 F ( D ( 1 1 — D L ( 1 1 .LT.VIGC  TO  5 
IF(DU(I1-D (I I.LT.V1G0  TO  6 
GC  TO  A 

5 J*J*1 
IB ( I >* J 

CO  TO  A 

6 J*J*1 
iC(l»«J 

A CONTINUE 

JX*J 
MN*JP+1 
11*0 
156  J*0 

OC  22I*1,MM 
11*1-1 
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IFd.EQ.lt  TH*TQ 
IFd.EQ.lt  GO  TO  77 
TM«GGdlt 

IFCGGdlt.LE.  BLdll)  GO  TO  22 
ID  (11 ) *1 
77  J*J-»1 
TMJI»TH 

22  CCNTINUE 

IFUX.EQ.JXN.AND.J.EQ.JINt  GO  TO  14 

J 1N=J 

JXN*JX 

DO  23  K-l.IP 
DT ( K t *D (K  1 
D(Kt*DlKt-*DEL 
J*0 

DC  24  I =1 »HM 
11*1-1 

IFd.EQ.l)  GO  TO  777 
IF1GGC  111  .LE.BLdltt  GO  TO  24 
777  J*J*1 

1 F ( J6  .ME  .0  1 GQ  TO  778 
DBJd,K»*CBJ(ll> 

776  DK(  J,K  l = (DBJd,Kl-TNUd/DEL 

24  CCNTINUE 

23  D(Kl*DT(Kt 
J6*l 

DO  25  1=1 ,J 

stm=o. 

DC  26  K*1 , I P 

26  SUM=SUM*CKd  ,Kl**2 

IF ( SUM .E Q .0 • t GO  TO  25 
5UM*D SORT (SUM  ) 

CO  27  K = 1 * I P 

27  CKd,K)=CKd,Kl/SUM 

25  CONTINUE 
JC=0 

330  IX=W*2 

1 F ( J.GT.l )G0  TO  775 
1FI JX.EG.OtGO  TO  150 
775  DO  350  1*1, LX 
340  DC  350  J=1 , B 
350  A ( I 1 *0 
J*0 

00  2 1*1, MM 
11*1-1 

IFd.EQ.lt  GO  TO  7 
IF  (GG  (1 1 1 »LE  . BL(  11 1 > GO  TO  2 
7 J*J*1 

L*0 

CO  3 K*1 , I P 

IFdCtKl.NE.OI  GO  TO  3 
L*L*1 

A ( J , L ) *DK  ( J ,K  1 
3 CCNTINUE 

DC  10  K*1,IP 
IFde(Kt.NE.O)  GO  TO  10 
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L*L+1 

A ( J t L ) =— DK l J , K ) 

10  CONTINUE 
N*L  + 1 

IF  ( J.EC.l  )G0  TO  2 
A ( J VN ) *CK (11) 
IF(KY.NE.0»A(J,NI*10. 

2 CONTINUE 

A(1,N>*1. 

J A*J+1 
OC  8 1*1, L 
J*J*1 

8 A ( J, I ) *1 >0 

W*J 

A (W+l «N)*1 *0 
B*J+N+1 
BB*B+1 
OC  9 I =JA » W 
9 A ( I »B I =1 .0 
L*W 

00  510  1*1, N 

510  A(U+l,n*FLOAT(Zl*A(W«l,Il 
M*J 
NA*6 

CALL  L INPRO (K2  1 
IFIK2.E0.0I  STOP 
LL*M+1 

OC  1460  1*1,  LL 
J=A ( I ,BB ) 

SSSi J)*Atl tBl 
1460  CENT 1NUE 
L *0 

DC  19  1*1,1P 
IFUCm.NE.OlGO  TO  12 
L *L+1 

5S<1 J=SSS(L) 

GC  TO  19 

12  SSUI*0 

19  CONTINUE 

OC  13  1*1,IP 

1 F ( IB ( I ) .N£ .0  I GO  TO  13 

L*L-*1 

ssm*ssm-sss(Li 

13  CCNTINUE 
14  SIIH*0 ■ 

DO  11  I-l.IP 

11  SUN*SUM«S5m«SS(I  I 
SUH*OSQRT (SUH  I 
IMSUH.EQ.O.IGO  TO  187 
OC  16  I*1,1P 

ssm*v*ss(  n/suH 

0(1 >*OT(l !*SS(I > 

16  CCNTINUE 

18  1K*0 

CALL  FINDIST) 
IFUP.EQ.OIGQ  T0157 
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J*0 

00  185  I*1,JP 
IF  CTC I I.LE.0.1  GO  tQ  185 
1F(ID< II .EQ.O)  BLU)*-1.5«DABS(GG(  I)  J 
IFdDm.EQ.Ol  J*1 
185  CONTINUE 

IFU.GT.OIGQ  TQ  15 
157  1FISB.GT.ST)  GO  TO  101 
187  V=.5*V 

CO  113  1*1, JP 
113  BL(Il*BL{Il/2. 

1FIV.LT.VT)  RETURN 
1FU.LT.AM1N)  RETURN 


15  DC  163  1*1, IP 
163  CIlbDTlII 
€0  TO  117 
101  JK=0 
RETURN 

150  SUMM*0.0 

DO  1 Jl*i,lP 

1 SimM=SUMM+DKll,Jll«»2 
SUMH=DSQRT (SUHM I 
IFISUHM-EC.Ol  JK  = 2 
IFISUMM.EQ.O.I  RETURN 
162  DC  151  1*1, IP 

D ( 1 1*DT ( I )-DK (1,1 1«»/SUMH 
IFIDtlUT.DLdll  DtlhDLU) 

151  IFIDm.GT.DUtIM  D(l)-DUU) 
SUMH*0. 

DC  IbZ  1*1, IP 

152  5UMM*SUMM*(D(II-DTCin*«2 
5UMM=DSQRTl SUMM 1 

DC  153  1*1, IP 

153  D(!)*DT(11+(D(1!-DT(1!)4V/SUHM 

IK*0 

CALL  F IND (ST) 

IF( 5ti.GT.ST IGO  TQ  101 
V*V/2. 

CO  112  1*1, JP 
112  BL(n*8L(l)/2. 

lF(V.LT.VT)  RETURN 
1F(V.GT.AH1N)G0  TO  162 
RETURN 
E 80 
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SUBROUTINE  SIZE 

IMPLICIT  REAL*8 <A-h,Q-Z> 

COMMON  D(10),P(100),DZ(10I,X(IO),G(10),A,AMIN,DL(10),DU(10), 
IV (55,6b  I (CK (10  I ,0,IP,JP,1A ( 10 ) ,KP ,KW,KL  , I K , 1 DUM (22  1 
TR*OBJ (01 
DC  1 1*1, IP 
X(I)*D(I) 

0(1) =0(11 +.00 1 

0Z(I )*(0BJ(0)-TN)/.001 

1 0(1) *X (I) 

DT  EMP  =0.0 
1TEMP*0 
DC  2 1*1,  IP 

1F(DABS(0Z(I)).LT.0TEMP)  GO  TO  2 
CTEMP*DABS(DZ( I 1 ) 

ITEMP*I 

2 CONTINUE 
A*OABS((0.01«TN)/DZ(1 TEMPI) 

IFU.LT.O  .005  JA=0-005 

SUM  =0 • 

DC  3 1*1, IP 

3 5UM*5UM+D(I)*»2 
SUM *D SORT (SUM  1 4 .001 
IF(A.LT.SUM)  A*SUM 
AMI N*A /100000  « 

RETURN 

E AO 
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SUBROUTINE  LINPRQIK21 

IMPLICIT  REAL*8IA-H,G-ZI 
INTEGER  Z,E,G,BB,W,B,R,C,H 

CCHMON  DI10J,FC1001,0Z  1 10 ) ,X 1 10 J ,T 1 10 ) , V, AM  I N.DL 1 1C  I ,DU 1 10 1 . 

1 A (55, 66  I,  CM  10) ,Q,IP,JP,IA (10 ) , KP ,K W,Kl , I K , I L ,M ,N ,B  ,L  ,88. E »G,Z ,H , 
2U  » IOUM  (10  1 ,JX 
K2*l 
NA*6 
M*M-1.0 
560  LL*M*2 

OC  580  K=2,LL 
570  AtK-l.N+G+K-l )*1 
580  A ( K-l , BB 1 *K+N*G-1 
600  1 F (G.N6.G  1 GO  TO  620 

610  IF(E.EQ.O)  GO  TO  780 

611  GC  TO  650 
620  KK*L+E*2 

LL*M*2 

OC  630  K*KK  »L L 

630  AIK-1,K+N-L-E-1)*-1  ’ , 

650  h*N4l 

660  0*0  . i • 

670  LL*N*G  -(. 

OC  760  J*1 ,LL 
680  S*0 

690  LL1*M-G-E*2 

KK1*M+1 

OC  700  I =LL1 , KK1 
700  S*S4A(I,JI 
720  AIW*1,JI*-S 

730  lF(AtU*l»Jl.CT.Qk  GO  TO  760 
790  QSA(W+1 , J 1 
750  C*J 
760  CCNTINUE 
761  S*0 

762  LL*M-G-E*2 
KK*M+1 

DC  763  J*LL,KK 
763  S*S*A(J,B1 
765  A (ta+l,B I*-S 
780  CONTINUE 
790  IF(G.EQ.O)  GO  TO  810 
LL*N+1 
KK=N+G 

810  IF(L.EQ.O)  GO  TO  830 
LL*N4G*1 
KK*N+G-*L 
IC*1 

830  lFIG+E.EQ.O)  GO  TO  2000 
831  Ll*N4G4L*l 
KK*B-1 
1C*1 

860  GO  TO  2000 

895  IFIQ.EQ.99999)  GO  TO  1230 

900  IF(Q.EQ.O)  GO  TO  1330 
910  GC  TO  1900 
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920  H*H+1 

930  Q*.1E39 

940  R*-l 

LL*H+1 

950  DC  1000  1 *1 «LL 

960  IF(All,C> .LE.O)  GO  TO  1000 

970  1FCA(I  ,B)/A(1,C>.GT.Q)  GO  TO  1000 

980  Q=A(I,B)/AII,C) 

990  R*  I 
1000  CONTINUE 

1010  IF  l FLOAT! R 1 .GE  .-.5  I GO  TO  1050 
1 C = 2 

1030  GO  TO  2000 
1050  P*A (R ,C  ) 

1060  A(R,BBl*C 
1070  DO  1080  J*1,B 
1080  A(R,J1=A(R,J1/P 

LL=W*1 

1100  DC  1180  1*1 » LL 

1110  lF(l.EQ.R)  GO  TO  1180 

1120  DC  1170  J*1,B 

1130  lFIJ.EQ.Cl  GO  TO  1170 

1140  A(ltJlsA(IiJ)*-A(RiJ)*A(IiC) 

1150  IF(0ABS(A(1,J1I.GT..1E-4I  GO  TO  1170 

1160  A ( 1 « JI =0 
1170  CONTINUE 
1180  CONTINUE 
LL=W+1 

1190  DO  1200  1*1,  LL 
1200  A ( I ,C 1 *0 
1210  CONTINUE 
1220  A ( R ,C  I *1 
1230  Q*0 

1240  Ll*N+G+L 

DC  1280  J*1 »LL 

1250  1F(AIW+1»JI.&T.Q)GC  TO  1280 

1260  Q*A(N+1,J1 

1270  C = J 

1280  CONTINUE 

1290  GG  TO  900 

1330  1F(W.EQ.M+1)  GO  TO  1360 

1340  H =W- 1 

1350  IF(A(W*2,B1.LT..1E-5I  GO  TO  1353 

K2*0 

1352  RETURN 

1353  LL*H+1 

DO  1358  1*1,  LL 

1354  IF(IDlNT(Atl,6Bll.LE.N«G«U  GO  TO  1358 

1355  DC  1356  J*1,B 

1356  A ( 1 , J ) *0 

1358  CONTINUE 

1359  GO  TO  1230 
136C  CONTINUE 

1400  IFIQ.EQ.OI  GO  TO  1420 

1420  LL*M«l 

1430  DC  1460  1*1 ,LL 
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1AA0  lF(IDINT(A(ltBB)).EQ.O)  GO  TO  1460 

1460  CONTINUE 

1470  IF(Q.NE.O)  GO  TO  920 

1520  KK*N«1 

LL«B-G-1 
XJJ*-Z«A(W*1,B  ) 

Lt*H-l 
I C*3 

1550  GC  TO  2000 
2000  L L*H-1 
Ll*W*l 

GC  TO(895«1C50  »999  I ( 1 Q 
RETURN 
END 


999 


- ■ 


OBJ 


r 


REAL  FUNCTIUN  OB J«8 I J 1 
IMPLICIT  REAL*B(A-H»Q-Z) 

CCMMQN  X(10 J,P(100»,DUH(3693) , IDUH138I 
101  OBJ>B-U 

IF  (B.EQ.G  . .OR.U.EQ.O.)  RETURN 
OBJ*lB-U)/DABS(Ul 
RETURN 
END 
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